### This script plots the analyses to compare different ways to incorporate the PDG
### into the Zonation analyses
### Libraries
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 3.5.2
library(dplyr)
## Warning: package 'dplyr' was built under R version 3.5.2
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(tidyr)
## Warning: package 'tidyr' was built under R version 3.5.2
### Read data
## Read output data generated by the script spatial_analyses_zonationVSproxiesdivgen_all_WT_final.R
## this shows the area per species and proxies of genetic diversity within it
## considering the Zonation solution for 20% of the country considering only
## species distribution models (no habitat, etc).
sols_summary_spp<-read.csv("../data/comparations_output/sols_summary_spp.txt")
# check data
head(sols_summary_spp)
## Solution Richness shannon simpson Area Prop_to_AreaSP
## 1 ENM_sp 64 3.216956 0.9341559 473887 1.0000000
## 2 Zonation_MDP 58 2.908383 0.9179912 127904 0.2699040
## 3 Zonation_MDP_ZV 59 2.929987 0.9210733 127082 0.2681694
## 4 Zonation_MDP_PDG 61 3.565616 0.9677778 142830 0.3014010
## 5 Zonation_MDP_vs_PDG 64 3.718850 0.9720398 124221 0.2621321
## 6 Zonation_MDP_PDG_ADMU 63 3.182794 0.9418937 129665 0.2736201
## mean.prop median.prop Dist.to.Optimun sp
## 1 NA NA 0 Capsicum annuum glabriusculum
## 2 0.3786268 0.2580767 NA Capsicum annuum glabriusculum
## 3 0.4083096 0.3092418 NA Capsicum annuum glabriusculum
## 4 0.5480085 0.4495986 NA Capsicum annuum glabriusculum
## 5 0.6525528 0.6421699 NA Capsicum annuum glabriusculum
## 6 0.4407066 0.3329768 NA Capsicum annuum glabriusculum
### Analyses and plots
## Check mean proportion of species area and mean proportion of proxies within it, per solution
group_by(sols_summary_spp, Solution) %>%
summarise(., mean.prop.sp.area = mean(Prop_to_AreaSP) , mean.prop.proxies = mean(mean.prop))
## # A tibble: 6 x 3
## Solution mean.prop.sp.area mean.prop.proxies
## <fct> <dbl> <dbl>
## 1 ENM_sp 1 NA
## 2 Zonation_MDP 0.476 0.543
## 3 Zonation_MDP_PDG 0.367 0.576
## 4 Zonation_MDP_PDG_ADMU 0.481 0.663
## 5 Zonation_MDP_vs_PDG 0.408 0.757
## 6 Zonation_MDP_ZV 0.429 0.555
### Plots
## Violin plots
# dist to optimun
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Solution, y=Dist.to.Optimun, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
stat_summary(fun.y=mean, geom="point", color="red") +
theme(axis.text.x = element_blank()) +
scale_y_continuous(name="Distance to Optimun")
## Warning: Removed 595 rows containing non-finite values (stat_ydensity).
## Warning: Removed 595 rows containing non-finite values (stat_summary).
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in min(diff(sort(x))): no non-missing arguments to min; returning
## Inf
## Warning in min(x): no non-missing arguments to min; returning Inf
## Warning in max(x): no non-missing arguments to max; returning -Inf
## Warning in stats::runif(length(x), -amount, amount): NAs produced
## Warning: Removed 595 rows containing missing values (geom_point).
# Simpson
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Solution, y=simpson, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
stat_summary(fun.y=mean, geom="point", color="red") +
theme(axis.text.x = element_blank()) +
scale_y_continuous(name="Simpson Index")
# Area sp conserved
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Solution, y=Prop_to_AreaSP, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
stat_summary(fun.y=mean, geom="point", color="red") +
theme(axis.text.x = element_blank()) +
scale_y_continuous(name="Proportion of distribution of the sp conserved")
# Mean prop area proxies by sp ponserved
plot_b <-filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Solution, y=mean.prop, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
stat_summary(fun.y=mean, geom="point", color="red") +
theme(axis.text.x = element_blank()) +
scale_y_continuous(name="Mean proportion of the area of each proxy conserved by sp")
plot_b
# Median prop area proxies by sp ponserved
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Solution, y=median.prop, color=Solution)) + geom_violin() + geom_jitter(size = 0.3) +
stat_summary(fun.y=mean, geom="point", color="red") +
theme(axis.text.x = element_blank()) +
scale_y_continuous(name="Median proportion of the area of each proxy conserved by sp")
# Median and mean together
filter(sols_summary_spp, Solution!="ENM_sp") %>%
gather(., SummaryStat, value, -Solution, -c(Richness:Prop_to_AreaSP), -Dist.to.Optimun, -sp) %>%
ggplot(., aes(x=interaction(SummaryStat, Solution), y=value, color=SummaryStat)) + geom_violin() + geom_jitter(size = 0.1) +
theme(axis.text.x = element_text(angle = 45, hjust = 1))
# Proportion of the distribution of the sp conserved vs MEAN proportion of each proxy conserved by sp
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Prop_to_AreaSP, y=mean.prop, color=Solution)) + geom_point() +
scale_x_continuous(name="Proportion of distribution of the sp conserved") +
scale_y_continuous(name="Mean proportion of the area of each proxy conserved by sp")
# Proportion of the distribution of the sp conserved vs MEADIAN proportion of each proxy conserved by sp
filter(sols_summary_spp, Solution!="ENM_sp") %>%
ggplot(., aes(x=Prop_to_AreaSP, y=median.prop, color=Solution)) + geom_point() +
scale_x_continuous(name="Proportion of distribution of the sp conserved") +
scale_y_continuous(name="Median proportion of the area of each proxy conserved by sp")
## with regression lines
plot_a <-filter(sols_summary_spp, Solution!="ENM_sp") %>%
# plot points
ggplot(., aes(x=Prop_to_AreaSP, y=mean.prop, color=Solution)) + geom_point(size=0.7) +
scale_x_continuous(name="Proportion of species distribution (%)") +
scale_y_continuous(name="Mean proportion of area of proxy represented by species distribution (%)",
breaks = seq(0, 1, by = 0.25), expand = c(0, 0)) +
# plot fitted curve
geom_smooth(method=loess, aes(fill=Solution))
plot_a
### Plot for paper figure:
plot_a <- plot_a +
# nicer background
theme_bw()+ theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
# nicer legend (fill and color are needed becasue we are plotting both points smooth and points)
scale_fill_discrete(name = "Scenarios",
labels = c("SDM (n=116)",
"SDM+LZ (n=143)", "SDM+PGD (n=218)",
"SDM*PGD (n=5004)", "SDM and PGD as ADMU (n=117)")) +
scale_color_discrete(name = "Scenarios",
labels = c("SDM (n=116)",
"SDM+LZ (n=143)", "SDM+PGD (n=218)",
"SDM*PGD (n=5004)", "SDM and PGD as ADMU (n=117)")) +
theme(legend.position = c(.99, 0.01), legend.justification = c("right", "bottom")) +
# title
labs(title="a)") +
# larger text
theme(
plot.title = element_text(size=14, face="bold"),
axis.title = element_text(size=14),
axis.text = element_text(size=13),
legend.text = element_text(size=13),
legend.title = element_text(size=13, face="bold"))
plot_a
plot_b <- plot_b +
# nicer background
theme_bw()+ theme(panel.border = element_blank(), panel.grid.major = element_blank(),
panel.grid.minor = element_blank(),
axis.line = element_line(colour = "black")) +
# no legend
theme(legend.position = "none") +
# nicar x names
scale_x_discrete(name="Scenarios",
labels=c("SDM",
"SDM+LZ", "SDM+PGD",
"SDM*PGD", "as ADMU")) +
# title
labs(title="b)") +
# larger text
theme(
plot.title = element_text(size=14, face="bold"),
axis.title = element_text(size=14),
axis.text = element_text(size=13))
plot_b
## Plot Together
# Multiple plot function taken from http://www.cookbook-r.com/Graphs/Multiple_graphs_on_one_page_(ggplot2)/
#
# ggplot objects can be passed in ..., or to plotlist (as a list of ggplot objects)
# - cols: Number of columns in layout
# - layout: A matrix specifying the layout. If present, 'cols' is ignored.
#
# If the layout is something like matrix(c(1,2,3,3), nrow=2, byrow=TRUE),
# then plot 1 will go in the upper left, 2 will go in the upper right, and
# 3 will go all the way across the bottom.
#
multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) {
library(grid)
# Make a list from the ... arguments and plotlist
plots <- c(list(...), plotlist)
numPlots = length(plots)
# If layout is NULL, then use 'cols' to determine layout
if (is.null(layout)) {
# Make the panel
# ncol: Number of columns of plots
# nrow: Number of rows needed, calculated from # of cols
layout <- matrix(seq(1, cols * ceiling(numPlots/cols)),
ncol = cols, nrow = ceiling(numPlots/cols))
}
if (numPlots==1) {
print(plots[[1]])
} else {
# Set up the page
grid.newpage()
pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout))))
# Make each plot, in the correct location
for (i in 1:numPlots) {
# Get the i,j matrix positions of the regions that contain this subplot
matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE))
print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row,
layout.pos.col = matchidx$col))
}
}
}
multiplot(plot_a, plot_b, cols=2)
